Introduction

Types de données

Il y a 3 principaux types de données spatiales

Données géostatistiques
Données de points (point pattern)
Données d’aires (areal data or lattice data

Données géostatistiques

Ce sont des données mesurées à des points fixes dans l’espace

Aussi connu sous le nom de géostatistiques, kriging.

Données récoltées à des stations météo

geoR, gstat, nlme, glmmTMB, R-INLA

Point pattern

Ce sont des données où c’est la localisation des points dans l’espace qui est aléatoire

L’ensemble des mentions d’occurrence pour la Phragmite commune (Phragmites autralis)

Upside-down sloths are so cute spatstat, R-INLA

Areal data

Ce sont des données qui sont agrégées au niveau d’entités surfaciques

Proportions d’oiseaux testant positifs au virus du Nil par unités administratives

Upside-down sloths are so cute spdep, R-INLA

Un exemple de cas

Il s’agit d’un jeu de données concernant le nombre d’espèces de plantes richesse (nb d’espèces) à différents parcelles sur l’île Las Palmas dans les îles Canaries. Cette île est caractérisée par de fortes variations en altitude en raison de la présence de deux volcans sur l’île. Ces données proviennent d’une étude par Irl et al. (2015) et peuvent être téléchargées librement à cette adresse:

https://doi.org/10.5061/dryad.7sk2g/1

Les exemples qui suivent sont fortement inspirés de Zuur et al. (2017).

La base de données légèrement modifiée peut être téléchargée directement à partir de la page github de l’atelier.

d <- read.csv("https://raw.githubusercontent.com/frousseu/introRINLA/master/canary_richness.csv")
head(d)

##   Plotname  easting northing richness altitude  TCI
## 1    AC001 221.1035 3190.789       21   0.7784 1.26
## 2    AC002 221.4035 3192.189       19   0.4473 1.14
## 3    AC003 221.8035 3192.689       25   0.2921 1.28
## 4    AC005 222.1035 3192.889       21   0.0610 1.48
## 5    AC006 222.0035 3192.689       23   0.1424 1.78
## 6    AC007 222.4035 3192.989       22   0.1211 1.86

On a donc un identifiant de parcelle, des coordonnées, la richesse en espèce et deux variables. L’altitude est en kilomètre et TCI est un index de complexité topographique.

En premier, utilisons les fonctions du package sp pour convertir le data.frame en objet spatial. Ensuite, nous téléchargeons l’île en question avec la fonction getData du package raster.

library(sp)
library(raster)
library(scales)

coordinates(d) <- ~easting + northing  # convert to spatial
proj4string(d) <- "+proj=utm +zone=28 +ellps=WGS84 +datum=WGS84 +units=km +no_defs"  # assign projection
ds <- spTransform(d, CRS("+init=epsg:4326"))  # transform to latlon

spa <- getData("GADM", country = "ESP", level = 2)  # download Spain as a shapefile
laspalmas <- disaggregate(spa)[225, ]  # subset Las Palmas

par(mar = c(3, 3, 0, 0))  # plot locations ans Las Palmas
plot(laspalmas, axes = TRUE)
plot(ds, add = TRUE, col = alpha("darkgreen", 0.5), pch = 16)  # alpha adds tranparency to points

On peut également visualiser les localisation avec une carte interactive à l’aide du package mapview.

library(mapview)
mapviewOptions(basemaps = c("Esri.WorldImagery", "Esri.WorldShadedRelief"))
mapview(ds, zcol = "richness")

Avant de faire un modèle, on visualise brièvement les données. Le @data permet d’extraire le data.frame (la table d’attributs) de l’objet d qui est maintenant un objet spatial, plus précisément un SpatialPointsDataFrame.

plot(d@data[, c("richness", "altitude", "TCI")])

Utilisons premièrement un simple GLM avec un effet quadratique sur l’altitude.

fit <- glm(richness ~ altitude + I(altitude^2) + TCI, data = d@data, family = poisson)
summary(fit)

## 
## Call:
## glm(formula = richness ~ altitude + I(altitude^2) + TCI, family = poisson, 
##     data = d@data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.8273  -1.3622  -0.2198   0.8298   6.8234  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.89824    0.04299  44.159   <2e-16 ***
## altitude       0.58252    0.06330   9.203   <2e-16 ***
## I(altitude^2) -0.58296    0.03322 -17.550   <2e-16 ***
## TCI            0.67431    0.02908  23.186   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4493.8  on 889  degrees of freedom
## Residual deviance: 2813.5  on 886  degrees of freedom
## AIC: 6655.6
## 
## Number of Fisher Scoring iterations: 5

On peut visualiser le modèle avec le package visreg. L’argument scale permet d’obtenir les prédictions ous l’échelle de la réponse (et non sous l’échelle log).

library(visreg)
par(mfrow = c(1, 2), mar = c(4, 4, 0, 1))
visreg(fit, "altitude", scale = "response")
visreg(fit, "TCI", scale = "response")

Maintenant, plusieurs des parcelles sont très près les unes des autres et on peut se demander si les nombres d’espèces observés sont plus similaires pour les parcelles près les unes des autres. En d’autres mots, reste-t-il de la variation inexpliquée par nos variables explicatives qui pourrait être structurée spatialement? Pour déterminer ceci, on peut faire un variogramme à l’aide des fonctions du package geoR.

Pour ce faire, on peut premièrement extraire les résidus (variation non-expliquée) de notre modèle et les cartographier afin de déterminer visuellement s’il y a des patrons suggérant la présence de corrélation spatiale. Si les petits résidus ou les grands résidus ont tendance à être ensemble, cela nous suggère qu’il y a potentiellement de la corrélation spatiale.

par(mar = c(3, 3, 0, -0))
plot(d, cex = rescale(resid(fit), to = c(0.2, 3)), col = gray(0.5, 0.5), pch = 16, 
    axes = TRUE)

Il semble y avoir un tel patron, mais cela n’est pas si facile à confirmer visuellement. On peut vérifier cela de façon plus formelle à l’aide d’un variogramme.

Intuitivement, un variogramme représente la variance des différences entre toutes les paires d’observations pour différentes classe de distances.

library(geoR)
v <- variog(coords = cbind(d$easting, d$northing), data = resid(fit))

## variog: computing omnidirectional variogram

plot(v, main = "Empirical Variogram for Species Richness", type = "b", xlab = "Distance (km)", 
    ylab = "Semivariance")

Si les observations près dans l’espace se ressemble davantage, on devrait s’attendre à ce que les valeurs de semivariance soient plus faibles pour les distances courtes. En fait, une courbe plate suggère qu’il n’y a pas de corrélation spatiale, alors qu’une courbe qui augmente (et qui atteint possiblement un plateau) suggère qu’il y a corrélation.

Nous avons pris les valeurs par défaut de la fonction variog et notamment le nombre de classes de distances qui est relativement faible. Augmentons le nombre de classes afin d’avoir un peu plus de précision sur ce qui se passe à faible distance.

v <- variog(coords = cbind(d$easting, d$northing), data = resid(fit), breaks = seq(0, 
    20, by = 0.5), max.dist = 20)

## variog: computing omnidirectional variogram

plot(v, main = "Empirical Variogram for Species Richness", type = "b", xlab = "Distance (km)", 
    ylab = "Semivariance")

On peut voir que la variance est beaucoup moins élevée à de faibles distances. On peut contrôler davantage le variogramme pour inspecter la forme de la relation.

v <- variog(coords = cbind(d$easting, d$northing), data = resid(fit), breaks = seq(0, 
    10, by = 0.5), max.dist = 10)

## variog: computing omnidirectional variogram

plot(v, main = "Empirical Variogram for Species Richness", type = "b", xlab = "Distance (km)", 
    ylab = "Semivariance")

Il semble donc y avoir une importante corrélation jusqu’à environ 2km, après quoi les points sont plutôt plats, indiquant qu’au delà de cette distance, la corrélation est faible ou inexistante.

On a donc une corrélation spatiale et il faut donc en tenir compte dans notre analyse! Pour cela, il faut utiliser des méthodes spéciales qui sont disponibles dans les packages geoR et gstat. Cependant, ces packages peuvent être utilisés avec des modèles gaussiens, mais pas avec les GLM (à part geoRglm).

INLA facilite grandement l’inclusion d’effets spatiaux notamment dans les GLMs ou GLMMs auxquels nous sommes habitués.

Variogrammes

Différentes valeurs de paramètres

Variogramme empirique

Un variogramme empirique est un variogramme généré à partir de données.

## variog: computing omnidirectional variogram

Variogramme théorique

Les variogrammes théoriques sont les modèles qui sont ajustés aux variogrammes empiriques.

library(gstat)
show.vgms(models = c("Exp", "Gau", "Sph", "Cir", "Mat"), as.groups = TRUE)

Différents types de variogrammes:

Exp: exponentiel
Gau: gaussien
Sph: sphérique
Cir: circulaire
Mat: Matérn

Les différents types de variogrammes sont décrits par des paramètres qui peuvent être ajustés.

show.vgms(kappa.range = c(0.5, 1, 2, 5), models = "Mat", nugget = c(0.1, 0.2, 
    0.3, 0.4), max = 10, as.groups = TRUE)

Un variogramme est souvent décrit par:

range (étendue de la corrélation spatiale)
sill (le plateau atteint en variance)
nugget (la variation due à d’autres facteurs)

v <- show.vgms(range = 2, nugget = 3, sill = 6, models = "Exp", max = 10, as.groups = TRUE, 
    ylim = c(0, 11), plot = FALSE)
plot(semivariance ~ distance, data = v[-1, ], type = "l", ylim = c(0, 10), col = "blue", 
    xaxs = "i", yaxs = "i", xlab = "distance", ylab = "semivariance")
abline(6, 0, lty = 3)
abline(9, 0, lty = 3)
abline(3, 0, lty = 3)
abline(v = 2, lty = 3)

INLA

INLA veut dire Integrated Nested Laplace Approximation. Pour une description “accessible” de la théorie derrière INLA, je vous suggère ce lien dans un blogue.

Approche bayésienne

Bayes’ rule

\[P(A|B) = \frac{P(B|A)P(A)}{P(B)}\]

Lorsque appliquée aux modèles

\[P(\theta|data) = \frac{P(data|\theta)P(\theta)}{P(data)}\]

\(data\): données \(\theta\): paramètres du modèle

Priors

LGM

SPDE

SPDE est pour Stochastic Partial Differential Equation

Cette approche se base entre autres sur une discrétisation du Gaussian Random Field (GRF) par un Gaussian Markov Random Field (GMRF)

Utilise une fonction de Matérn pour capturer la dépendance entre les localisations.

Tiré de Zuur, Ieno et Saveliev (2017), Volume I

library(INLA)
m <- inla(richness ~ altitude + I(altitude^2) + TCI, data = d@data, family = "poisson")
summary(m)

## 
## Call:
## c("inla(formula = richness ~ altitude + I(altitude^2) + TCI, family = \"poisson\", ",  "    data = d@data)")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          2.3899          0.7240          0.3430          3.4568 
## 
## Fixed effects:
##                  mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept)    1.8980 0.0430     1.8138   1.8980     1.9826  1.8978   0
## altitude       0.5828 0.0633     0.4589   0.5826     0.7074  0.5823   0
## I(altitude^2) -0.5831 0.0332    -0.6487  -0.5830    -0.5182 -0.5827   0
## TCI            0.6744 0.0291     0.6169   0.6745     0.7311  0.6747   0
## 
## The model has no random effects
## 
## The model has no hyperparameters
## 
## Expected number of effective parameters(std dev): 4.078(0.00)
## Number of equivalent replicates : 218.25 
## 
## Marginal log-Likelihood:  -3348.80

m$summary.fixed[, c(1:3, 5)]

##                     mean         sd 0.025quant 0.975quant
## (Intercept)    1.8980388 0.04300921  1.8137909  1.9825871
## altitude       0.5827786 0.06333174  0.4588786  0.7074154
## I(altitude^2) -0.5830944 0.03323408 -0.6486835 -0.5182087
## TCI            0.6743690 0.02910015  0.6168778  0.7311290

cbind(summary(fit)$coef[, 1:2], confint(fit))

##                 Estimate Std. Error      2.5 %     97.5 %
## (Intercept)    1.8982431 0.04298664  1.8142114  1.9827229
## altitude       0.5825224 0.06329847  0.4589618  0.7070942
## I(altitude^2) -0.5829604 0.03321702 -0.6484524 -0.5182393
## TCI            0.6743125 0.02908278  0.6169035  0.7309101

Un exemple de cas

Cas

Étapes

L’estimation d’un modèle avec un effet spatial avec INLA requiert de passer par plusieurs étapes.

Mesh: permet de créer une triangulation qui discrétise le GF
Projector Matrix: permet de relier les localisations à la grille
SPDE: permet détablir le lien le SPDE et la fonction Matérn
Spatial Field: permet de créer le champs spatial
Stack: permet d’intégrer les différents éléments
Formula: formulation du modèle
Modèle: permet de faire tourner le modèle

Voici un résumé graphique de ces différentes étapes avec les fonction associées.

Mesh

La première étape consiste à créer une grille (mesh) qui va être utilisée pour approximer le champs Gaussien. En premier, il faut obtenir les localisations.

locs <- coordinates(d)
head(locs)

##    easting northing
## 1 221.1035 3190.789
## 2 221.4035 3192.189
## 3 221.8035 3192.689
## 4 222.1035 3192.889
## 5 222.0035 3192.689
## 6 222.4035 3192.989

Par la suite, on utilise la fonction inla.mesh.2d pour créer la grille. Optionnellement, on peut utiliser une grille non-convexe qui respecte davantage le contour de nos observations.

library(INLA)
hull <- inla.nonconvex.hull(locs, convex = -0.05)
mesh <- inla.mesh.2d(loc = locs, offset = c(1, 5), max.edge = c(1, 3), cutoff = 1, 
    boundary = hull)

par(mar = c(0, 0, 6, 0))
plot(mesh, asp = 1)
points(locs, col = alpha("black", 0.6), pch = 16, cex = 0.7)

offset
Spécifie l’étendue de l’extension de la grille au delà des observations et l’extension d’une zone tampon au-delà de cette grille. Cette dernière permet de réduire les effets de bordure lors des estimations. L’étendue de cette zone tampon doit être au moins aussi grande que l’étendue de la corrélation spatiale.

max.edge
Permet de spécifier les dimensions maximales des triangles de la grille à l’intérieur et à l’extérieur dans la zone tampon. Plus ces triangles sont petits, plus les approximations sont précises, mais plus les calculs sont longs. En général, il n’est pas nécessaire que les triangles dans la zone tampon soient aussi petits qu’à l’intérieur de la grille.

cutoff
Par défaut, chaque observation sera utilisée comme un coin d’un triangle (vertex). Pour éviter la création de trop nombreux petits triangles, on spécifie une valeur à cutoff au delà de laquelle des points voisins seront ignorés lors de la création des triangles.

boundary
Cet argument permet d’utilier l’étendue non-convexe autour des observations.

Il faut éviter d’avoir des triangles avec des aigles trop aigus. Autrement, les estimations sont de moins bonnes qualités.

Projector Matrix

Ceci permet de relier les localisations à la grille et d’établir une pondération qui permet d’estimer les valeurs du champs spatial pour chaque localisation. Les localisations situées sur les coins (vertex) seront estimées à partir des valeurs de la grille et les localisations à l’intérieur du triangle seront estimées à partir d’une moyenne pondérée calculée en utilisant les trois coins du triangle dans lequel elles se trouvent.

A <- inla.spde.make.A(mesh, locs)

Dans cet exemple, l’observation en rouge est située à l’intérieur d’un tiangle et les valeurs représentent la pondération qui sera utilisée sur les valeurs de chaque coin pour estimer la valeur de ce point.

par(mar = c(0, 0, 5, 0))
plot(mesh, asp = 1, xlim = c(223, 226), ylim = c(3172, 3176))
points(locs, pch = 16, cex = 2)
i <- 889  # ligne contenant la localisation en rouge
points(locs[i, , drop = FALSE], col = "red", pch = 16, cex = 2)
w <- which(A[i, ] > 0)  # pondérations associées au point
points(mesh$loc[w, 1:2], label = 1:nrow(mesh$loc), pch = 1, cex = 4, lwd = 2, 
    col = "red")
text(mesh$loc[w, 1:2], label = round(A[i, w], 2), adj = c(-0.6, 0.5), font = 2, 
    col = "red")

SPDE

Ceci permet de définir les éléments du SPDE et les éléments associés aux caractéristiques du champs spatial.

spde <- inla.spde2.pcmatern(mesh, alpha = 2, prior.range = c(2, 0.5), prior.sigma = c(5, 
    0.01))

L’argument alpha est pour spécifier une des paramètres de la fonction de Matérn. Ce paramètre doit être fixé entre 0 et 2. Nous devons également définir les priors du champs spatial. L’approche par défaut est d’utiliser des quantiles pour définir les priors sur l’étendue (range) et l’écart-type (standard deviation, sd) associés au champs spatial. Cette façon de faire se base sur l’approche des Penalized-complexity priors développée par Simpson et al. (2017). En résumé, cette approche permet de pénaliser l’étendue de la corrélation spatiale vers l’infini (ce qui réduit la complexité du phénomène) et la variance vers 0 (ce qui réduit également la complexité du phénomène). En pratique, le prior pour l’étendue \(r\) de la corrélation est spécifié avec \(\alpha\) et \(r_{0}\):

\[P(r<r_{0})= \alpha\]

Où \(\alpha\) représente la probabilité que \(r\) soit inférieur à \(r_{0}\). Dans notre cas, prior.range=c(2,0.5) veut dire que nous croyons qu’il y a 50% de chances que le l’étendue de la corrélation spatiale soit supérieure à 2km et donc qu’il y a également 50% des chances q’elle soit inférieure à 2km. Le prior pour la variance du champs spatial, \(\sigma\), est spécifié avec \(\alpha\) et \(\sigma_{0}\):

\[P(\sigma>\sigma_{0})= \alpha\]

Où \(\alpha\) représente la probabilité que \(\sigma\) soit inférieur à \(\sigma_{0}\). Dans notre cas, cela indique qu’on croit qu’il y a 1% des chances que la variance spatiale en richesse soit supérieure à 5. Ici, il faut nous rappeler que nous travaillons sous l’échelle log et que \(e^5\approx148\) espèces.

Spatial Field

Ceci permet de définir un index qui sera utile pour spécifier ou récupérer les éléments associés au champs spatial.

spatial.index <- inla.spde.make.index(name = "spatial", n.spde = spde$n.spde)

Stack

Le stack est une façon compliquée de fournir les données, les variables et les effets à INLA. Ceci n’est pas essentiel pour des modèles simples, mais ça le devient lorsque les modèles sont plus compliqués ou lorsque l’on veut par exemple générer des prédictions à partir de notre modèle.

d$altitude2 <- d$altitude^2
v <- c("TCI", "altitude", "altitude2")
X <- data.frame(Intercept = 1, d@data[, v])
stack <- inla.stack(data = list(richness = d$richness), A = list(A, 1), effects = list(spatial = spatial.index, 
    as.list(X)), tag = "est")

Formula

On créé la formule décrivant le modèle souhaité. Notez que l’intercept est créé à la mitaine afin de.

model <- richness ~ -1 + Intercept + altitude + altitude2 + TCI + f(spatial, 
    model = spde)

Modèle

Finalement, on fait tourner le modèle en appelant la fonction inla et en fournissant les différents arguments. La spécification de compute=TRUE fait en sorte que les distributions postérieures seront calculées pour toutes les observations. La spécification de config=TRUE permettra plus loin de faire des simulations et de générer des valeurs aléatoires à partir de notre modèle.

m <- inla(model, data = inla.stack.data(stack), control.predictor = list(A = inla.stack.A(stack)), 
    family = "poisson")
summary(m)

## 
## Call:
## c("inla(formula = model, family = \"poisson\", data = inla.stack.data(stack), ",  "    control.predictor = list(A = inla.stack.A(stack)))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          6.8725         23.1234          0.5559         30.5519 
## 
## Fixed effects:
##              mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## Intercept  1.9972 0.1018     1.7965   1.9973     2.1967  1.9977   0
## altitude   0.6313 0.1694     0.2989   0.6311     0.9646  0.6307   0
## altitude2 -0.5258 0.0804    -0.6837  -0.5259    -0.3678 -0.5260   0
## TCI        0.4624 0.0511     0.3617   0.4625     0.5623  0.4628   0
## 
## Random effects:
## Name   Model
##  spatial   SPDE2 model 
## 
## Model hyperparameters:
##                     mean     sd 0.025quant 0.5quant 0.975quant   mode
## Range for spatial 2.1257 0.3607     1.4873   2.1031      2.902 2.0639
## Stdev for spatial 0.4277 0.0272     0.3772   0.4267      0.484 0.4245
## 
## Expected number of effective parameters(std dev): 216.06(7.583)
## Number of equivalent replicates : 4.119 
## 
## Marginal log-Likelihood:  -2894.36

Visualisations

Le modèle généré (m) est un objet complexe contenant énormément d’information où il n’est pas toujours facile de s’y retrouver. Heureusement, comme nous venons de le voir la fonction summary peut être utilisée pour un sommaire rapide. Si on veut extraire plus d’informations sur notre modèle, la première étape est probablement d’inspecter les différents éléments de la liste formant m.

names(m)

##  [1] "names.fixed"                 "summary.fixed"              
##  [3] "marginals.fixed"             "summary.lincomb"            
##  [5] "marginals.lincomb"           "size.lincomb"               
##  [7] "summary.lincomb.derived"     "marginals.lincomb.derived"  
##  [9] "size.lincomb.derived"        "mlik"                       
## [11] "cpo"                         "po"                         
## [13] "waic"                        "model.random"               
## [15] "summary.random"              "marginals.random"           
## [17] "size.random"                 "summary.linear.predictor"   
## [19] "marginals.linear.predictor"  "summary.fitted.values"      
## [21] "marginals.fitted.values"     "size.linear.predictor"      
## [23] "summary.hyperpar"            "marginals.hyperpar"         
## [25] "internal.summary.hyperpar"   "internal.marginals.hyperpar"
## [27] "offset.linear.predictor"     "model.spde2.blc"            
## [29] "summary.spde2.blc"           "marginals.spde2.blc"        
## [31] "size.spde2.blc"              "model.spde3.blc"            
## [33] "summary.spde3.blc"           "marginals.spde3.blc"        
## [35] "size.spde3.blc"              "logfile"                    
## [37] "misc"                        "dic"                        
## [39] "mode"                        "neffp"                      
## [41] "joint.hyper"                 "nhyper"                     
## [43] "version"                     "Q"                          
## [45] "graph"                       "ok"                         
## [47] "cpu.used"                    "all.hyper"                  
## [49] ".args"                       "call"                       
## [51] "model.matrix"

Distributions postérieures

Les premiers éléments qui devraient nous intéresser sont les distributions postérieures des coefficients associés à nos différentes variables. On peut extraire les quantiles de celles-cià partir du sommaire des effets fixes.

m$summary.fixed

##                 mean         sd 0.025quant   0.5quant 0.975quant
## Intercept  1.9971786 0.10184005  1.7965115  1.9973494  2.1967014
## altitude   0.6313093 0.16943925  0.2989201  0.6310948  0.9645592
## altitude2 -0.5258307 0.08042505 -0.6837271 -0.5258781 -0.3678250
## TCI        0.4624021 0.05107928  0.3617380  0.4625302  0.5622610
##                 mode          kld
## Intercept  1.9976856 7.381960e-07
## altitude   0.6306687 6.299124e-07
## altitude2 -0.5259636 3.673678e-07
## TCI        0.4627912 2.249087e-06

On peut également représenter graphiquement ces distributions à l’aide des distributions marginales complètes. Comme on peut le voir, tous les coefficients sont relativement loin de zéros.

par(mfrow = c(2, 2), mar = c(4, 4, 1, 1))
invisible(lapply(names(m$marginals.fixed), function(i) {
    p <- m$marginals.fixed[[i]]
    plot(p[, 1], p[, 2], type = "l", xlab = i, ylab = "density")
    abline(v = 0, lty = 3)
}))

Paramètres spatiaux

Le sommaire des paramètres associés au champs spatial s’accèdent par le sommaire des “hyperparamètres” (hyperparameters). Avec INLA, les paramètres fixes représentent les paramètres de la régression et tous les paramètres de variances ou associés au champ spatial sont considérés comme des hyperparamètres.

m$summary.hyperpar

##                        mean         sd 0.025quant  0.5quant 0.975quant
## Range for spatial 2.1257011 0.36066113  1.4872934 2.1031495  2.9020523
## Stdev for spatial 0.4277356 0.02722437  0.3772038 0.4267148  0.4839872
##                        mode
## Range for spatial 2.0639373
## Stdev for spatial 0.4245023

On peut également illustrer les distributions postérieures de ces paramètres.

par(mfrow = c(1, 2), mar = c(4, 4, 1, 1))
invisible(lapply(names(m$marginals.hyperpar), function(i) {
    p <- m$marginals.hyperpar[[i]]
    plot(p[, 1], p[, 2], type = "l", xlab = i, ylab = "density")
}))

Champ spatial

Le champs spatial peut être visualiser de différentes façons. La façon la plus simple est d’utiliser les fonction de INLA pour projeter les valeurs du champs sur une grille. Ceci peut être fait de cette façon.

library(fields)
library(viridisLite)

# https://haakonbakka.bitbucket.io/btopic108.html#92_plotting_the_spatial_mean_field


xlim <- range(d$easting)
ylim <- range(d$northing)

# - Can project from the mesh onto a 300x300 plotting grid
proj <- inla.mesh.projector(mesh, xlim = xlim, ylim = ylim, dims = c(300, 300))

# - Do the projection
mfield <- inla.mesh.project(projector = proj, field = m$summary.random[["spatial"]][["mean"]])
sdfield <- inla.mesh.project(projector = proj, field = m$summary.random[["spatial"]][["sd"]])


par(mfrow = c(1, 2), mar = c(3, 3, 2, 3))

image.plot(list(x = proj$x, y = proj$y, z = mfield), col = viridis(100), asp = 1)
axis(1)
axis(2)

image.plot(list(x = proj$x, y = proj$y, z = sdfield), col = viridis(100), asp = 1)
axis(1)
axis(2)

Prédictions

Les prédictions dans INLA sont générées en soumettant des observations pour lesquelles la variable réponse est NA (dans notre cas, richness=NA). Si on utilise l’argument compute=TRUE, INLA se chargera de calculer des distributions postérieures pour l’ensemble des valeurs avec NA. Cela peut nous permettre par exemple de soumettre des valeurs d’altitude pour étudier son effet sur la richesse.

On génère 50 valeurs d’altitude entre les valeurs minimales et maximales.

n <- 50
x <- seq(min(d$altitude), max(d$altitude), length.out = n)
newX <- data.frame(Intercept = 1, TCI = mean(d$TCI), altitude = x, altitude2 = x^2)

Pour l’exercice, on va assumer que l’effet spatial est constant en prenant un endroit quelconque sur l’île. On répète cette localité autant de fois que le nombre de valeurs d’altitude. Idéalement, il serait plus logique de ne pas tenir compte de l’effet spatial dans les prédictions, mais cela complexifie passablement le code nécessaire.

newlocs <- matrix(c(224, 3175), ncol = 2)[rep(1, n), ]

On associe la localisation fictive à la grille.

A.pred <- inla.spde.make.A(mesh = mesh, loc = newlocs)

On construit une stack pour fournir les valeurs à prédires et on la joint à celle générée plus haut qui contient les données en tant que tel. On spécifie tag=pred ce qui va nous permettre d’identifier les lignes contenant les prédictions.

stack.pred <- inla.stack(data = list(richness = NA), A = list(A.pred, 1), effects = list(spatial = spatial.index, 
    as.list(newX)), tag = "pred")
stack.full <- inla.stack(stack, stack.pred)

Avec ce tag=pred, on construit un index qui nous permettra de récupérer les lignes contenant les prédictions.

index.pred <- inla.stack.index(stack.full, tag = "pred")$data

Finalement, on fait tourner notre modèle en prenant soin de spécifier compute=TRUE pour que INLA calcule les distributions postérieures. On spécifie également link=1 pour s’assurer que INLA retransforme les valeurs sous l’échelle du nombre d’espèces et non sous l’échelle log.

m <- inla(model, data = inla.stack.data(stack.full), control.predictor = list(A = inla.stack.A(stack.full), 
    compute = TRUE, link = 1), family = "poisson")

On récupère les valeurs prédites à partir du sommaire des valeurs prédites (summary.fitted.values) en utilisant l’index créé plus haut.

p <- m$summary.fitted.values[index.pred, c("0.025quant", "mean", "0.975quant")]

plot(x, p[, "mean"], ylim = range(unlist(p)), type = "l", xlab = "Altitude en km", 
    ylab = "Richesse")
lines(x, p[, "0.025quant"], lty = 3)
lines(x, p[, "0.975quant"], lty = 3)

Exercice 1

Liens utiles

Livres en ligne

Livres

Articles importants

Forum R-INLA

Select the menu on the left to expand / collapse table of contents (TOC) entries. Press button below to collapse all TOC except the top level headings.

If you only want to collapse level 3 headings press this button.

Introduction aux modèles spatiaux avec R-INLA

François Rousseu

3 mai 2018

Introduction

Types de données

Données géostatistiques

Point pattern

Areal data

Un exemple de cas

Variogrammes

Variogramme empirique

Variogramme théorique

INLA

Approche bayésienne

Priors

LGM

SPDE

Un exemple de cas

Cas

Étapes

Mesh

Projector Matrix

SPDE

Spatial Field

Stack

Formula

Modèle

Visualisations

Distributions postérieures

Paramètres spatiaux

Champ spatial

Prédictions

Exercice 1

Liens utiles

Livres en ligne

Livres

Articles importants

Forum R-INLA